Jak wszyscy wiedzą, z mediów społecznościowych można się o każdym dowiedzieć wszystkiego. Powszechny dostęp do wielu z tych serwisów został w ciągu ostatniego roku mocno ograniczony (np. Facebook) albo też obwarowany zdobyciem serii kodów dostępu etc. Na potrzeby dzisiejszych zajęć zajmiemy się Twitterem; dostęp do niego jest możliwy dzięki bibliotece rtweet, niestety, prawdopodobnie nie uda się jej uruchomić na kompuetrach w Pracowni... Twitter wymaga aż czterech różnych kodów, aby można było korzystać z API - te poniżej zostały stworzone specjalnie na potrzeby zajęć TEXT i po około tygodniu zostaną zresetowane.
# PRZYKŁAD 10.1
library(rtweet)
create_token(
app = "TEXTclass",
consumer_key = "mjBZhJWhlcYYkjTNQX1zC9z1t",
consumer_secret = "J8KjngzwuimbpanNzGuWJOhx6flNID5D9wyhguCpxI5HMzjWTt",
access_token = "414634084-T9H7IJR7nX0k2wYuJE8SWydgyaPduV4L0CweL3M1",
access_secret = "EKQBhBnAf6Ew9Q9adeY1bn7YJEykKdsYpUbTUuyknp4ZW")
## <Token>
## <oauth_endpoint>
## request: https://api.twitter.com/oauth/request_token
## authorize: https://api.twitter.com/oauth/authenticate
## access: https://api.twitter.com/oauth/access_token
## <oauth_app> TEXTclass
## key: mjBZhJWhlcYYkjTNQX1zC9z1t
## secret: <hidden>
## <credentials> oauth_token, oauth_token_secret
## ---
Po tym uwierzytelnaniu można już popbrać dane - poniżej zbierzemy 10000 tweetów w konta CNN
# PRZYKŁAD 10.2
tw <- get_timeline("cnn", n = 10000)
tw
## # A tibble: 3,234 x 88
## user_id status_id created_at screen_name text source
## <chr> <chr> <dttm> <chr> <chr> <chr>
## 1 759251 10845368. 2019-01-13 19:45:05 CNN "GOP. Socia.
## 2 759251 10845330. 2019-01-13 19:30:05 CNN "DNA. Socia.
## 3 759251 10845292. 2019-01-13 19:15:07 CNN Peop. Socia.
## 4 759251 10845255. 2019-01-13 19:00:07 CNN Amer. Socia.
## 5 759251 10845217. 2019-01-13 18:45:09 CNN Puer. Socia.
## 6 759251 10845179. 2019-01-13 18:30:05 CNN "Rep. Socia.
## 7 759251 10845141. 2019-01-13 18:15:07 CNN Pres. Socia.
## 8 759251 10845104. 2019-01-13 18:00:06 CNN Walt. Socia.
## 9 759251 10845066. 2019-01-13 17:45:05 CNN Miss. Socia.
## 10 759251 10845034. 2019-01-13 17:32:23 CNN A ma. Socia.
## # . with 3,224 more rows, and 82 more variables: display_text_width <dbl>,
## # reply_to_status_id <chr>, reply_to_user_id <chr>,
## # reply_to_screen_name <chr>, is_quote <lgl>, is_retweet <lgl>,
## # favorite_count <int>, retweet_count <int>, hashtags <list>,
## # symbols <list>, urls_url <list>, urls_t.co <list>,
## # urls_expanded_url <list>, media_url <list>, media_t.co <list>,
## # media_expanded_url <list>, media_type <list>, ext_media_url <list>,
## # ext_media_t.co <list>, ext_media_expanded_url <list>,
## # ext_media_type <chr>, mentions_user_id <list>,
## # mentions_screen_name <list>, lang <chr>, quoted_status_id <chr>,
## # quoted_text <chr>, quoted_created_at <dttm>, quoted_source <chr>,
## # quoted_favorite_count <int>, quoted_retweet_count <int>,
## # quoted_user_id <chr>, quoted_screen_name <chr>, quoted_name <chr>,
## # quoted_followers_count <int>, quoted_friends_count <int>,
## # quoted_statuses_count <int>, quoted_location <chr>,
## # quoted_description <chr>, quoted_verified <lgl>,
## # retweet_status_id <chr>, retweet_text <chr>,
## # retweet_created_at <dttm>, retweet_source <chr>,
## # retweet_favorite_count <int>, retweet_retweet_count <int>,
## # retweet_user_id <chr>, retweet_screen_name <chr>, retweet_name <chr>,
## # retweet_followers_count <int>, retweet_friends_count <int>,
## # retweet_statuses_count <int>, retweet_location <chr>,
## # retweet_description <chr>, retweet_verified <lgl>, place_url <chr>,
## # place_name <chr>, place_full_name <chr>, place_type <chr>,
## # country <chr>, country_code <chr>, geo_coords <list>,
## # coords_coords <list>, bbox_coords <list>, status_url <chr>,
## # name <chr>, location <chr>, description <chr>, url <chr>,
## # protected <lgl>, followers_count <int>, friends_count <int>,
## # listed_count <int>, statuses_count <int>, favourites_count <int>,
## # account_created_at <dttm>, verified <lgl>, profile_url <chr>,
## # profile_expanded_url <chr>, account_lang <chr>,
## # profile_banner_url <chr>, profile_background_url <chr>,
## # profile_image_url <chr>
Pakiet rtweet ma wbudowaną funkcję ts_plot(), która korzysta z pakietu ggplot2 i do której, podobnie jak do ggplota, można dodawać kolejne warstwy. Argumentem funkcji ts_plot jest rozdzielczość z jaką chcemy wykreślić aktywnośc (czyli liczbę tweetów) w danych.
# PRZYKŁAD 10.3
ts_plot(tw, "5 hour", color = "gray") + geom_point(color = "blue")
Biblioteka oferuje całą masę różnych opcji (polecam prezentację oraz dokumentację). My zrobimy tylko jeden prosty wykres, ale wykorzystamy do niego wcześniej zdobyte metody - poliyczymy walencję w kolejnych tweetach, uzywając klasyfikatora słownikowego, opartego na zbiorze Warriner et al. z Laboratorium 6.
# PRZYKŁAD 10.4
library(tidytext)
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
emo <- as_tibble(read.csv("http://www.fizyka.pw.edu.pl/~julas/TEXT/lab/Ratings_Warriner_et_al.csv", stringsAsFactors = F))
emo <- emo %>%
select(word = Word, valence = V.Mean.Sum, arousal = A.Mean.Sum)
tw.data <- tw %>% transmute(id = row_number(), text = text) %>% unnest_tokens(word, text)
tw.sent <- tw.data %>%
inner_join(emo) %>%
group_by(id) %>%
summarise(valence = mean(valence))
## Joining, by = "word"
tw$valence[tw.sent$id] <- tw.sent$valence
## Warning: Unknown or uninitialised column: 'valence'.
ggplot(tw, aes(x = created_at, y = valence)) +
geom_line(alpha=0.25, color = "red") +
geom_smooth(method = "loess", span = 0.01)
## Warning: Removed 15 rows containing non-finite values (stat_smooth).
Google jest prawdziwym potentatem i w zasadzie monopolistą jeśli chodzi o szereg usług -- możemy albo z nich nie korzystać albo też zgadzać się na różne praktyki. Nie da się jednak ukryć, że poniższe usługi, do których dedykowane są konkretne API są dość przydatne:
Nas interesuje w tym momencie szczególnie ta ostatnia pozycja. R oferuje kilka pakietów do obsługi Google Maps - my skorzystamy z ggmap. Jak sama nazwa wskazuje, pakiet korzysta z ggplot2 i jest z nim kompatybilny pod względem idei użycia warstw. Na począwek wygenerujemy prostą mapę scentrowaną na środku Warszawy. Do obsługi map potrzebny jest nam klucz API - poniższy jest wygenerowany na użytek zajęć i zostanie zresetowany po dwóch tygodniach. Największym problemem jet to, że pakiet ggmap trzeba zainstalować z odrębnego źródła, a nie z ogólnego repozytorium i dodatkowo zajmuje to sporo czasu, bo niektóre pakiety są instalowane od nowa ze źródła.
# PRZYKŁAD 10.5
if(!requireNamespace("devtools")) install.packages("devtools")
devtools::install_github("dkahle/ggmap", ref = "tidyup")
# PRZYKŁAD 10.6
library(ggmap)
## Google Maps API Terms of Service: https://cloud.google.com/maps-platform/terms/.
## Please cite ggmap if you use it: see citation("ggmap") for details.
register_google("AIzaSyCVaYyyqIPta_LCss-P-P73k2fqAZPp7Ps")
map <- get_googlemap(center = c(lon = 21.00, lat = 52.25), zoom = 11, maptype = "road")
## Source : https://maps.googleapis.com/maps/api/staticmap?center=52.25,21&zoom=11&size=640x640&scale=2&maptype=roadmap&key=xxx-P-P73k2fqAZPp7Ps
ggmap(map)
Na mapę naniesiemy teraz pozycje przystanków sieci transportu miejskiego - zbiór wejściowy zosta pobrany ze serwera ZTM.
# PRZYKŁAD 10.7
stops <- readLines("http://www.if.pw.edu.pl/~julas/TEXT/lab/stops.txt")
head(stops)
## [1] " 1001 KIJOWSKA, -- WARSZAWA"
## [2] " *PR 9"
## [3] " 100101 2 Ul./Pl.: TARGOWA, Kier.: AL.ZIELENIECKA, Y= 52.248670 X= 21.044260 "
## [4] " L 6 - stały: 125 135 138 166 509 517 "
## [5] " L 1 - na żądanie: N21 "
## [6] " 100102 2 Ul./Pl.: TARGOWA, Kier.: ZĄBKOWSKA, Y= 52.249020 X= 21.044540 "
Na oko zbiór ma dość nieuporządkowaną postać, na szczęście okazuje się, że poszczególne dane, które chcielibyśmy na nasze potrzeby otrzymać (numer przystanku oraz współrzędne) łatwo poddają się ekstrakcji za pomocą prostych wyrażeń regularnych
# PRZYKŁAD 10.8
get.coor <- function(line) {
p <- regexpr("[0-9]{6}", line)
y <- regexpr("Y= 52.[0-9]{5,6}", line)
x <- regexpr("X= 2[0-1].[0-9]{5,6}", line)
p <- substr(line, p[1], p[1] + attr(p, "match.length") - 1)
y <- substr(line, y[1] + 3, y[1] + attr(y, "match.length") - 1)
x <- substr(line, x[1] + 3, x[1] + attr(x, "match.length") - 1)
data.frame(stop_id = as.numeric(p), long = as.numeric(y), lat = as.numeric(x))
}
stops <- get.coor(stops[grep("Y= ", stops)])
head(stops)
## stop_id long lat
## 1 100101 52.24867 21.04426
## 2 100102 52.24902 21.04454
## 3 100103 52.24900 21.04398
## 4 100104 52.24990 21.04173
## 5 100105 52.25035 21.04386
## 6 100106 52.25001 21.04371
Korzystając z powyższej ramki danych, możemy teraz nanieść poszczególne punkty na mapę, dodając je jako kolejną warstwę ggplota.
# PRZYKŁAD 10.9
ggmap(map) + geom_point(data = stops, aes(x = lat, y = long), size=0.5, shape = 19)
## Warning: Removed 1286 rows containing missing values (geom_point).
Na bazie takich danych można się pokusić o np. naniesienie warstwy z interpolowaną gestością przystanków (warstwa stat_density2d). Przy okazji dokonamy również drobnej zmiany podstawowej mapy, ustawiając język na polski oraz kolor na czarno-biały.
# PRZYKŁAD 10.10
map <- get_googlemap(center = c(lon = 21.00, lat = 52.25), zoom = 11, maptype = "road", color = "bw", language = "PL")
## Source : https://maps.googleapis.com/maps/api/staticmap?center=52.25,21&zoom=11&size=640x640&scale=2&maptype=roadmap&language=PL&key=xxx-P-P73k2fqAZPp7Ps
ggmap(map) +
geom_point(data = stops, aes(x = lat, y = long), size=0.5, shape = 19) +
stat_density2d(data = stops, aes(x = lat, y = long, fill = stat(level)), geom = "polygon", alpha = 0.05, bins = 30) + scale_fill_gradient(low = "blue", high = "red")
## Warning: Removed 1286 rows containing non-finite values (stat_density2d).
## Warning: Removed 1286 rows containing missing values (geom_point).
Pozycje przystanków są na pewno interesujące, ale zwykle więcej uwagi poświęca się danym dynamicznym, zmieniającym się w czasie. Takie możliwości oferuje np. serwis Otwarte Dane pod egidą Urzędu Miasta Warszawy. Rejestrując tam konto otrzymujemy klucz API do różnych danych - nas interesuje serwis udostępniający w czasie rzeczywistym pozycje tramajów w Warszawie.
Do opracowania wyników zapytań potrzebujemy biblioteki jsonlite, która obsługuje format JSON zwracany przez API
# PRZYKŁAD 10.11
library(jsonlite)
##
## Attaching package: 'jsonlite'
## The following object is masked from 'package:rtweet':
##
## flatten
um.waw.api <- "41f78414-bc00-4099-b607-42bb80a9d9d5"
um.waw.url <- "https://api.um.warszawa.pl/api/action/wsstore_get/?id=c7238cfe-8b1f-4c38-bb4a-de386db7e776&apikey="
url.api <- paste(um.waw.url, um.waw.api, sep = "")
trams <- fromJSON(url.api)
head(trams$result)
## Status FirstLine Lon Lines Time Lat
## 1 RUNNING 28 20.92525 28 2019-01-13T20:48:56 52.26143
## 2 RUNNING 6 20.94394 6 2019-01-13T20:48:57 52.29683
## 3 RUNNING 6 20.92996 6 2019-01-13T20:48:59 52.29200
## 4 RUNNING 27 20.97193 27 2019-01-13T20:48:59 52.27275
## 5 RUNNING 28 20.90473 28 2019-01-13T20:48:58 52.23906
## 6 RUNNING 35 20.93497 35 2019-01-13T20:48:57 52.27289
## LowFloor Brigade
## 1 FALSE 3
## 2 FALSE 6
## 3 FALSE 8
## 4 FALSE 3
## 5 FALSE 6
## 6 FALSE 7
Teraz swobodnie wykorzystujemy naszą mapę Google, aby nanieść na nią pozycje tramwajów.
# PRZYKŁAD 10.12
map <- get_googlemap(center = c(lon = 21.00, lat = 52.25), zoom = 12, maptype = "road", color="bw", language = "PL")
## Source : https://maps.googleapis.com/maps/api/staticmap?center=52.25,21&zoom=12&size=640x640&scale=2&maptype=roadmap&language=PL&key=xxx-P-P73k2fqAZPp7Ps
ggmap(map) +
geom_point(data = trams$result, aes(x = Lon, y = Lat))
## Warning: Removed 18 rows containing missing values (geom_point).
Mapę można trochę "uatrakcyjnić", nakładając na punkty numer linii oraz dodatkowo kolorując według tych numerów.
# PRZYKŁAD 10.13
trams$result$FirstLine <- as.numeric(gsub(" ", "", trams$result$FirstLine))
## Warning: NAs introduced by coercion
ggmap(map) +
geom_point(data = trams$result, aes(x = Lon, y = Lat, color = as.factor(FirstLine)), size = 4) +
geom_text(data = trams$result, aes(x = Lon, y = Lat, label = FirstLine), size = 3) +
scale_color_discrete(guide = FALSE)
## Warning: Removed 18 rows containing missing values (geom_point).
## Warning: Removed 19 rows containing missing values (geom_text).
Jak wiadomo, zanieczyszczenie powietrza (czyli popularny smog) stanowi coraz większy problem w Polsce - jest to wypadkowa natężonego ruchu pojazdów oraz ogrzewania domów za pomocą kiespkiej jakości paliw, przy czym ten drugi czynnik jest raczej decydujący. W rezultacie takiej sytuacji, pojawiają się coraz to nowe przedsiębiorstwa, oferujące dostęp do metod moniotorowania poziomu zanieczyszczenia. Jedną z takich form jest Airly, która udostępnia mapę oraz apki, umożliwiające ocene stanu powietrza w danej chwili, jak również obejrzenie 24-godzinnej historii takich pomiarów oraz prognoze na kolejną dobę. Airly dostarcza również API: po zarejestrowaniu się, otrzymujemy klucz dający umożliwiający wykonanie maks. 50 zapytań na minutę oraz do 1000 na dobę.
Naszym pierwszym krokiem będzie zdobycie danych 45 czujników z obszaru Warszawy.
# PRZYKŁAD 10.14
airly.api = "0oRd9wGrqtnhN91H1UaAhPwKjRhGrmxy"
airly.url1 <- "https://airapi.airly.eu/v2/installations/nearest?lat=52.25&lng=21.00&maxDistanceKM=30&maxResults=45&apikey="
query <- paste(airly.url1, airly.api, sep = "")
czujniki <- fromJSON(query)
head(czujniki)
## id location.latitude location.longitude address.country address.city
## 1 2191 52.25574 20.99337 Poland Warszawa
## 2 1074 52.24233 20.99543 Poland Warszawa
## 3 2350 52.25445 20.98454 Poland Warszawa
## 4 3190 52.24577 21.01822 Poland Warszawa
## 5 2525 52.23689 20.99439 Poland Warszawa
## 6 729 52.26425 20.98856 Poland Warszawa
## address.street address.number address.displayAddress1
## 1 Inflancka 4b Warszawa
## 2 aleja "Solidarności" 113 Warszawa
## 3 aleja Jana Pawła II 80 Warszawa
## 4 Sowia 4 Warszawa
## 5 Krochmalna 1 Warszawa
## 6 aleja Wojska Polskiego 20 Warszawa
## address.displayAddress2 elevation airly sponsor.id sponsor.name
## 1 Inflancka 4A 101.07 TRUE 22 Aviva
## 2 aleja Solidarności 113 114.64 TRUE 49 eurobank
## 3 Jana Pawła II 80 108.31 TRUE 22 Aviva
## 4 Sowia 88.35 TRUE 22 Aviva
## 5 Krochmalna 114.79 TRUE 10 WarszawaOddycha
## 6 aleja Wojska Polskiego 98.09 TRUE 10 WarszawaOddycha
## sponsor.description
## 1 Airly Sensor's sponsor
## 2 Airly Sensor's sponsor
## 3 Airly Sensor's sponsor
## 4 Airly Sensor's sponsor
## 5 Airly Sensor is part of action
## 6 Airly Sensor is part of action
## sponsor.logo
## 1 https://cdn.airly.eu/logo/Aviva_1538146740542_399306786.jpg
## 2 https://cdn.airly.eu/logo/eurobank.jpg
## 3 https://cdn.airly.eu/logo/Aviva_1538146740542_399306786.jpg
## 4 https://cdn.airly.eu/logo/Aviva_1538146740542_399306786.jpg
## 5 https://cdn.airly.eu/logo/WarszawaOddycha.jpg
## 6 https://cdn.airly.eu/logo/WarszawaOddycha.jpg
## sponsor.link
## 1 https://wiemczymoddycham.pl/
## 2 <NA>
## 3 https://wiemczymoddycham.pl/
## 4 https://wiemczymoddycham.pl/
## 5 <NA>
## 6 <NA>
Jak widać, każdy z czujników ma unikalny numer (id) oraz podaną lokalizajcę. Będziemy teraz chcieli odpytać każdy z czujników - naszym celem jest stan powietrza w danej chwili (reprezentowany przez tzn indeks ACQ), ale przy okazji otrzymamy też całą masę innych danych, które wykorzystamy później. W tym momencie stworzymy sobie ramkę danych, zawierającą wartości indeksu oraz pozycje czujników, a następnie korzystając z ggmap naniesiemy te dane na mapę.
# PRZYKŁAD 10.15
get.air.mes <- function(id, api.key) {
query <- sprintf("https://airapi.airly.eu/v2/measurements/installation?installationId=%d&apikey=%s", id, api.key)
print(id)
#Sys.sleep(1)
fromJSON(query)
}
pomiary <- lapply(czujniki$id, get.air.mes, api.key = airly.api)
## [1] 2191
## [1] 1074
## [1] 2350
## [1] 3190
## [1] 2525
## [1] 729
## [1] 337
## [1] 658
## [1] 2724
## [1] 6615
## [1] 6532
## [1] 6864
## [1] 956
## [1] 6961
## [1] 26
## [1] 3410
## [1] 3507
## [1] 2325
## [1] 6959
## [1] 2679
## [1] 2330
## [1] 7448
## [1] 817
## [1] 528
## [1] 6630
## [1] 2720
## [1] 3328
## [1] 27
## [1] 2143
## [1] 812
## [1] 2205
## [1] 3511
## [1] 6963
## [1] 2067
## [1] 1071
## [1] 650
## [1] 1091
## [1] 1066
## [1] 2129
## [1] 6954
## [1] 2214
## [1] 3465
## [1] 732
## [1] 1016
## [1] 6996
indexes <- sapply(1:length(pomiary), function(i) pomiary[[i]]$current$indexes$value)
data <- cbind(czujniki$location, indexes)
head(data)
## latitude longitude indexes
## 1 52.25574 20.99337 35.66
## 2 52.24233 20.99543 35.81
## 3 52.25445 20.98454 32.97
## 4 52.24577 21.01822 32.50
## 5 52.23689 20.99439 37.44
## 6 52.26425 20.98856 NA
ggmap(map) + geom_point(data = data, aes(x = longitude, y = latitude, color = indexes), size = 4) +
scale_color_gradient2("AQI", low = "green", high = "violet", mid = "red", midpoint = 85)
Jak wspomiano, nie są to jedyne dane, które przechowujemy w zmiennej pomiary. Dla każdego czujnika mamy dostęp do hisorii oraz do prognozy: wykorzystajmy tę pierwszą informację do policzenia korelacji pomiędzy pomiarami poszczególnych czujników.
# PRZYKŁAD 10.16
library(corrplot)
historia <- sapply(1:length(pomiary), function(i) sapply(pomiary[[i]]$history$indexes, function(x) x$value))
C <- cor(historia)
corrplot(C, type = "upper", na.label = "-", diag = F)
Jak widać, korelacje są raczej wysokie - tylko jeden czujnik (nr 11) jest wyraźnie inny. Dysponując takimi danymi można się pokusić o sprawdzenie jak korelacja zachowuje się w funkcji odległości. Potrzebujemy do tego funkcji wyznaczającej odegłości pomiędzy punktami określonami przez współrzędne geograficzne (wykorzystamy bibliotekę geosphere). Poza tym użyjemy bibliotek fields i Hmisc do binowania oraz prezentowania słupków niepewności.
# PRZYKŁAD 10.17
library(geosphere, quietly = T, warn.conflicts = F)
library(fields, quietly = T, warn.conflicts = F)
## Spam version 2.1-1 (2017-07-02) is loaded.
## Type 'help( Spam)' or 'demo( spam)' for a short introduction
## and overview of this package.
## Help for individual functions is also obtained by adding the
## suffix '.spam' to the function name, e.g. 'help( chol.spam)'.
##
## Attaching package: 'spam'
## The following objects are masked from 'package:base':
##
## backsolve, forwardsolve
library(Hmisc, quietly = T, warn.conflicts = F)
D <- distm(data[c(2, 1)], fun = distHaversine)/1000
d <- stats.bin(D[upper.tri(D)], C[upper.tri(C)])
x <- d$centers
y <- d$stats[2,]
yerr <- d$stats[3] / sqrt(d$stats[1,])
par(mfrow = c(1,2))
plot(D[upper.tri(D)], C[upper.tri(C)], pch = 19, cex = 0.5, col = "gray", tcl = 0.25, xlab = "odelglosc [km]", ylab = "korelacja")
errbar(x, y, y + yerr, y - yerr, add = T, col = "blue", errbar.col = "blue")
errbar(x, y, y + yerr, y - yerr, col = "blue", errbar.col = "blue", tcl = 0.25, xlab = "odelglosc [km]", ylab = "korelacja")
par(mfrow = c(1,1))